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Abstract. We report simulations of the inspiral and merger of binary neutron stars 
performed with WhiskyTHC, the first of a new generation of numerical relativity codes 
employing higher than second-order methods for both the spacetime and the hydro- 
dynamic evolution. We find that the use of higher-order schemes improves substan¬ 
tially the quality of the gravitational waveforms extracted from the simulations when 
compared to those computed using traditional second-order schemes. The reduced de¬ 
phasing and the faster convergence rate allow us to estimate the phase evolution of the 
gravitational waves emitted, as well as the magnitude of finite-resolution effects, with¬ 
out the need of phase- or time-alignments or rescalings of the waves, as sometimes 
done in other works. Furthermore, by using an additional unpublished simulation at 
very high resolution, we confirm the robustness of our high convergence order of 3.2. 


1. Introduction 


The inspiral and merger of binary neutron stars (BNSs) is one of the most promising 
sources of gravitational waves (GWs) fo r future ground-based laser-inte rferometer de¬ 
tectors such as LIGO, Virgo or KAGRA (ISathvaprakash & Schutzii2009 ). Because they 
can travel almost unscattered through matter, GWs carry valuable information from the 
deep core of the neutron stars (NSs) concerning the equation of state (EOS) of matter 
at supra-nuclear densities. Unfortunately, they are also extremely hard to detect, so that 
their identification and analysis requires the availability of analytical or semi-analytical 
GW templates. In turn, the validation and tuning of these models must be done by 
matching them with the predictions of fully non-linear numerical relativity (NR) calcu¬ 
lations, which represent the only means to describe accurately the late in spiral of BNS 


(Baiotti et al. 201(|; 


Bernuzzi et al.ll2014 


Baiotti et al. 20 


Bernuzzi et al 


to describe accurately tne late mspirai or bix: 
Bernuzzi et al.ll2012al : lHotokezaka etaD 12013 


20141) . 


While very high-quality NR waveforms of binary black hole me rgers are available 


e.g. (lAvlott et al J 120091 : Mroue & et al.ll2013l : iHinder & et al.ll2013h . BNS simulations 


have been pla gued by low convergence order and relat ively large phase uncertainties 
(Sf If ~ 1%) feaiotti et al.ll2009l : lBemuzzi et al.ll2012bh . Furthermore, since NSs have 


smaller masses, the merger part of the waveform is out of the frequency band for the 
next generation GW detectors, so that EOS-related effects will have to be most probably 
extracted from the inspiral signal using a large number of events. In particular, EOS- 
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induced effects will be encoded in the phase evolution of the GW signal during the 
inspiral dDamour et al.ll2012h . As a result, the measure of EOS-induced effects requires 


very accurate general-relativistic predictions of the inspiral signal. Even though accu- 


costs 

Baiotti et al. 2010; Baiotti et al.l 201 

Ill: Bernuzzi et al. 2012aj; Hotokezaka et al. 

2013; 

Bernuzzi et al. 20 1J; Bernuzzi et al. 

20l4), their analysis is complicated by the 


low convergence order of the methods employed. In particular, the analysis often re¬ 
quires the use of a time rescaling or alignm ent of the GWs from different resolutions 
dBaiotti et al.ll201ll : lHotokezaka et al.ll2013h . which is hard to justify mathematically. 

Here we show that, by using high-order numerical methods, it is indeed possible 
to obtain waveforms for the late-inspiral of a BNS system with a quality that is al¬ 
most comparable with the one obtained for binary black holes, with clean, higher than 
second-order convergence in both the phase and the amplitude. In particular we high¬ 
light and extend to higher resolution some of the results we reported in dRadice et al 
2014a). 


2. Numerical methods 


The results presented here have been obtained with our new high-order, high-resolution 
shock-capturing (HRSC), finite-differencing code: WhiskyTHC dRadice et al.ll2014al Jbh. 
which represents the extension to general relativity of the special-relativistic THC code 
dRadice & Rezzollal l2012h . Wh i s k y T HC solves the equa t ions of general-rela tivistic 
hydrodynamics in conservation form dBanvuls et al.lll997l : iRezzolla & Zanottih using 
a finite-difference scheme that employs flux reconstruction in local-characteristic vari¬ 
ables using the MP 5 scheme, formally fifth-orde r in space dSuresh & Huvnhl[l997h [see 
Radice & Rezzollal d2012h : iRadice et aD d2014bh for details]. 

The s pacetime evolution makes use of the BSSNOK formulation of the Einstein 
equations (iNakamura et al Jl 1 987l : IShibata & Nakamuralfl995l : lBaumgarte & Shapirol[l999b 
and it is performed using fourth-order accurate finite-differenc e sc heme_provided by the 


Me la chi an and is part of t he Einstein Toolkit dLoffler et al.ll2012l : lBrown et al. 


20091 : ISchnetter et ai1l2004) . To ensure the non-linear stability of the scheme we add 
a fifth-order Kreiss-Oliger type artificial dissipation to the spacetime variables only. 
Finally, the coupling between the hydrodynamic and the spacetime solvers is done us¬ 
ing the method of lines and a fourth-order Runge-Kutta time integrator. The resulting 
scheme is formally fourth-order in space and time, except at the boundaries between 
different refinement levels, where our method is only second-order in time. This should 
have only marginal effects on our results given that our finest grid covers both NSs. 


3. Binary setup 


The initial data is computed in the conformally flat approximation using the Lorene 
pseudo-spectral code ( Gourgoulhon et al.1 12001 ) and describes two irrotational, equal- 
mass NSs in quasi-circular orbit. Its main properties are summarized in Table [U and 
we note that it is computed using a polytropic EOS (p - Kp r ) with K = 123.56 (in 
units where G - M© - c - 1) and T = 2, while th e ev olution unperformed using the 
ideal-gas EOS (p = (T - 1 )pe) with the same T ( Rezzolla & Zanotti ). 
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Table 1. Summary of the considered BNS model. We report the total baryonic 
mass Mb, the ADM mass M, the initial separation r, the initial orbital frequency 
/orb, the gravitational mass of each star at infinite separation, Moo, the compactness, 
C = Moo/tfco, where Roc is th e area l radius of the star when isolated and the tidal 
Love number, a: 2 , e.g. iHinderer et al. ( 20101) . 


M b [M e ] 

M [M q \ 

r [km] 

/orb [HZ] 

Moo [Mq] C K2 

3.8017 

3.45366 

60 

208.431 

1.7428 0.18002 0.05 


The runs are performed on a grid covering 0 < x, z < 750 km, -750 km < y < 
750 km, where we assume reflection symmetry across the (x,y) plane and n symmetry 
across the (y, z) plane. The grid employs six fixed refinement levels, with the finest 
one covering both stars. We consider four different resolutions, labelled as L, M, H 
and VH , having, in the finest refinement level, a grid spacing of h ^ 370,295,215 and 
147 meters, r espec tively. The results of simulations L, M, H were already presented in 
Radice etldl d2014al lbh. while those of the simulation VH are presented here for the first 
time. 


4. Binary dynamics 

The dynamics of the binary is summarized in Figure [U where we show the density on 
the equatorial plane for the M run (the others are qualitatively very similar). The initial 
distance between the two NS centers is 60 km. They complete about 7 orbits, while 
inspiraling because of the loss of orbital angular momentum to GWs, before entering 
into contact at time t ^ 24.6 ms (from the beginning of the simulation). Finally, they 
quickly merge and a black hole is formed soon after. 


5. Gravitational waves 


To quantify the accuracy of our code for GW astronomy, we consider the phase evolu¬ 
tion of the dominant t - 2, m - 2 component of the curvature GW, the (complex) Weyl 
scalar ¥ 4 , as extracted at r ^ 665 km. We compute the GW phase / after decomposing 
the curvature as ¥4 = Ae -1< ^. Figure [2] shows an analysis of the residual of the phase 
between the different resolutions as a function of the retarded time u. On the left panel, 
we show both the absolute de-phasing between successive resolutions and the residuals, 
between the H and M and the VH and H resolutions, scaled assuming a convergence 
order of 3.2. On the right panel, we plot the instantaneous convergence order as mea¬ 
sured from three out of the four resolutions (separately the first three L, M, H and the 
last three M, H , VH). 

Figure [2] demonstrates that higher than second-order convergence can be achieved 
for numerical relativity simulations of BNS, without the need to perform any artificial 
manipulation of the waveforms. We find an order of convergence ~ 3.2 (also con¬ 
firmed by the very good overlap between the rescaled de-phasing), which is somewhat 
smaller than the formal order of four of our scheme. However, this is to be expected 
because HRSC methods typically reach their nomina l conver gence order only at very 
high resolutions (IShull 1 9971 : iRadice & Rezzollall2012h : see also lZlochower et all 
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Figure 1. Rest-mass density on the equatorial plane for the medium resolution 
run M. The panels show the initial data (upper right), the late stages of the inspiral 
(top middle and left), the approximate time of contact (bottom left), merger (bottom 
center) and black-hole formation (bottom right). 




Figure 2. Accumulate de-phasing (left panel) and estimated order of convergence 
(right panel) for the € = 2, m = 2 mode of the Weyl scalar ¥4 as extracted at 
r = 450 M q . The de-phasing between high ( H ) and medium (M) and very high (VH) 
and high are also rescaled assuming an order of convergence of 3.2. The instanta¬ 
neous order of convergence is estimated separately from the first three resolutions 
<9({L, Af, H}) and from the last three <9({M, H, VH}). 
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Figure 3. Relative phase (left panel) and amplitude (right panel) differences for 
the € = 2, m = 2 mode of the Weyl scalar ¥4 between the Richardson extrapolated 
phase and each of the resolutions. The extrapolated phase is computed assuming a 
convergence order of 3.2 and using the first three resolutions (L, M and H). 


for a discussion of other possible sources of errors. Finally, as also observed with other 
codes dBernuzzi et al.ll2012bh . our solution shows a loss of convergence (with apparent 
super-convergence) after u > 24.6 ms. This is roughly the time when the two NSs 
enter into contact (see Figure Q]). At this time the de-phasing between the H and VH 
resolution is - 0.26 rad. 

As a measure of the error with respect to the exact solution we use the first three 
resolutions (L, M and H ), Richardson extrapolate the GW to infinite resolution assum¬ 
ing a convergence order of 3.2 and compute the phase and amplitude differences of 
each run with respect to the extrapolated data. The results are shown in Figure [3j Ob¬ 
viously the data for u > 24.6 ms (which is about 13.5 GW cycles) has to be taken with 
a grain of salt, given that convergence is lost after contact. Before that, we find that the 
H run has a de-phasing ^ 0.35 rad (- 0.4%) at time u - 24.6 ms, while the highest 
resolution VH run has a de-phasing as small as ^ 0.13 rad. This corresponds to a phase 
error of less than 0.15%, which would be challenging to achieve with standard second- 
order numerical-relativity codes. The relative errors of the amplitude are somewhat 
larger, but also on a few percent level for the H and VH resolutions, which is more than 
adequate giving that the amplitude is not as critical as the phase for GW astronomy. 
Overall, this shows that the Richardson extrapolated waveform is consistent with the 
VH data (which lies between the extrapolated waveform and the H data). 


6. Conclusions 

We presented a set of four, high-resolution, high-order, simulations of BNS inspiral in 
general relativity, the highest of which has not been published before. Our analysis 
focused on the accuracy of the gravitational wave signal extracted from the simulations 
and, in particular, on its phase evolution, which represents a crucial quantity for both 
detection and parameters extraction. 

We showed that, with the use of higher-order methods, it is possible to achieve 
clean convergence and small phase errors, without the need to perform any alignment 
or rescaling of the GWs. The Richardson-extrapolated waveforms appear robust as 
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the resolution increases, giving support to their accuracy. This opens the possibility 
of obtaining high-quality waveforms with reliable error estimates that could be used to 
verify and calibrate analytical and phenomenological phasing models to be used in GW 
astronomy. 
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